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Abstract 

In this paper, we consider the problem of steering a family of independent, structurally identical, finite-dimensional stochastic 
linear systems with variation in system parameters between initial and target states of interest by using an open-loop control 
function. Our exploration of this class of control problems, which falls under the rising subject of ensemble control, is motivated 
by pulse design problems in quantum control. Here we extend the concept of ensemble control to stochastic systems with 
additive diffusion and jump processes, which we model using Brownian motion and Poisson counters, respectively, and consider 
optimal steering problems. We derive a Fredholm integral equation that is used to solve for the optimal control, which 
minimizes both the mean square error (MSE) and the error in the mean of the terminal state. In addition, we present several 
example control problems for which optimal solutions are computed by numerically approximating the singular system of the 
associated Fredholm operator. We use Monte Carlo simulations to illustrate the performance of the resulting controls. Our 
work has immediate practical applications to the control of dynamical systems with additive noise and parameter dispersion, 
and also makes an important contribution to stochastic control theory. 
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1 Introduction 

The behavior of physical, chemical, and biological sys- 
tems can exhibit significant sensitivity to uncertainty 
or variation in system parameters. This factor arises in 
practical control problems in many areas of science and 
engineering when there is uncertainty in the parame- 
ters of a single control system, or when a collection of 
structurally similar systems with variation in common 
parameters must be steered using a common control sig- 
nal. The former case has been approached by employing 
H x feedback con trol such as for quadratic stabilization 
of linear systems (|Xie et al.l . ll992f) . where parameter un- 
certainty is modeled with unknown time-dependent vari- 
ations in the system dynamics. When feedback is un- 
available, however, analysis of these cases has given rise 
to the subject of ensemble control, which is motivated by 
the practical application of control theory to the fields of 
Nuclear Magnetic Resonance (NMR) spectrosco py and 
imaging (MRI) (P. 120061 ILTand Khaneial . l2006f) . 

Rapidly developing technologies based on quantum 
theory require the manipulation of large quantum en- 
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sembles, e.g., a spin ensemble on the order of Avo- 
godro's number (6 x 10 23 ), whose individual states 
cannot be measured and whose dynamics are subject 
to dispersion in natural frequencies. Such applications 
require open-loop controls whose effect is insensitive to 
dispersion in system parameters as well as to inhomo- 
genei ties in the appli e d rad i o-freq uency (RF) control 
field (iLi and Khaneial I2006L 120091 ). A typical class of 
control problems significant to applications of NMR 
and MRI requires the design of RF pulses that drive a 
collection of spin systems with identical dynamics but 
parameter values unknown up to a given range between 
initial and target states, such that control performance is 
immune to physical para meter variation (jKobzar et all 
120051: iPaulvet all . Il99ll) . A similar scenario appears m 
the field of neuroscience, where population level fre- 
quency control of neural oscillators requires open-loop 
stimuli that minimize signal power while achieving the 
desired response across a c ollection of neurons wit h vari- 
ation in natura l dyna mics (jBrown et all 120041 : . 120101 
IZlotnik and Lit 1201 if ) . Although initially motivated by 
the need to control large collections of similar systems, 
the paradigm of ensemble control can be readily used 
to solve any open-loop control problem in which the 
system response must be insensitive to uncertainty in 
model parameters. 



Many practical control systems are subject to random 



dynamic disturbances, which are best treated from a 
probabilistic perspective by including a stochastic term 
in the system model. We focus on the influence on en- 
semble control design of additive noise processes that 
perturb the system dynamics. Because parameter vari- 
ation can greatly affect the behavior of stochastic sys- 
tems, especially when feedback is not used to attenuate 
disturbances, an analysis of ensemble control problems 
that include stochastic terms is crucial and meaningful. 

It is important to clarify the distinction between pa- 
rameter uncertainty and stochastic effects with respect 
to system dynamics. Variation, or dispersion, in system 
parameters among members of a class of dynamical sys- 
tems gives rise to differences in the actuation of indi- 
vidual members, each of which has fixed parameter val- 
ues. In practical applications, factors such as imprecise 
models or measurements and natural inhomogeneities 
result in uncertainty in parameters among a collection 
of systems, which is typically quantifiable within a given 
range. In contrast, the state trajectory of a physical sys- 
tem often varies from a nominal path due to dynamic 
effects caused by the environment or inherent to the sys- 
tem itself, which may change randomly and depend on 
time. These perturbations can be modeled as stochas- 
tic disturbances, whose nature may also depend on a 
parameter, so that the system state is described by a 
probability density function. For a given parameter set, 
such stochastic effects can be attenuated by minimizing 
statistical objectives such as the expectation and vari- 
ance of the system state distribution, as in the linear 
quadratic Gaussian (LQG) problem. Our contribution is 
an alternative method for designing open-loop controls 
that are insensitive to parameter dispersion while opti- 
mizing these statistical objectives. 

This paper is organized as follows. In the next section, 
we will provide a brief review of ensemble control for 
finite-dimensional time- varying deterministic linear sys- 
tems, and introduce the definition of stochastic ensemble 
controllability. In Section 3, we present our main result, 
which is a derivation of the conditions that the optimal 
ensemble control must satisfy in order to minimize the 
MSE at the terminal state. It is shown that the same 
control also minimizes the error in the mean, and co- 
incides with the optimal control for the corresponding 
deterministic case as well. We provide a derivation for 
the case of Brownian diffusion, and indicate the neces- 
sary modifications for considering a Poisson type jump 
process. In Section 4, we discuss the numerical approx- 
imation of solutions to the Fredholm integral equation 
of the first kind that characterizes the optimal control 
by discretizing the kernel in time and parameter space, 
so that the singular system can be approximated by the 
singular value decomposition (SVD) of a matrix opera- 
tor. Finally in Section 5, several examples are given in 
order to illustrate the application and performance of 
our method. 



2 Preliminaries 



In this section, we review the basic ideas of ensemble 
control and introduce the notion of stochastic ensem- 
ble control. Consider a parameterized family of finite- 
dimensional time- varying linear control systems 

X(t, (3) = A(t, /3)X(t, /3) + B(t, /3)u(t), 

X e M C R", (3 eK cR d , ueU C K m , (I) 

where A(t,P) £ R" x " and B(t,(3) G W ixm have ele- 
ments that are real and Li functions, respectively, 
defined on a compact set D = [0, T] x K C IR 2 , and are 
denoted A G L^ n (D) and B G L% xm (D). The ensem- 
ble controllability conditions for the system (1) depend 
on the existence of an open-loop control u ; [0, T] — > U 
that can steer the instantaneous state of the ensem- 
ble X : D — >• M between any points of interest. We 
consider the system (1) in a Hilbert space setting. Let 
Ht = L™[0, T] be the set of m-tuples, whose elements 
are vector- valued square- integrable functions defined on 
< t < T, with an inner product defined by 

(g,h) T = [ T g\t)h(t)dt, (2) 
Jo 

where ' denotes the transpose. Similarly, let Hk = 
LVi{K) be equipped with an inner product 

(p,q) K = I J{fi)qW)M0), ( 3 ) 

J K 

where fi is the Lebesgue measure, and K is the com- 
pact domain of the parameter j3. With well-defined ad- 
dition and scalar multiplication, Ht and Hk are sepa- 
rable Hilbert spaces, where 1 1 • | |t and \ \-\\k denote their 
respective induced norms. We now restate th e defin ition 
of deterministic ensemble controllability ([La , 1201 it) . 

Definition 1. We say that the family (1) is ensemble 
controllable on the function space Hk if for all e > 0, and 
all Xq, Xp G Hk , there exist T > and an open- loop 
piecewise-continuous control u G Ht, such that starting 
from X(0, f3) = X (p), the final state X(T, /3) = X T (f3) 
satisfies \\Xt — Xf\\k < e - 

In other words, the system is ensemble controllable if 
it is possible to drive the system from Xq to with 
respect to all j3 G K, where the acceptable range of 
T G (0, oo ) may depend on e, K, and U. Necessary and 
sufficient conditions have been provided for the ensemble 
controllability of finite-dimensional time- var ying linea r 
systems by applying integral operator theory (|Lil . l201ll) . 
Given the initial state X(0,f3) — X {(3) of the system 
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(1), the variation of constants formula gives rise to 

X{T,P) = $(T,O,0)X O (P) 

+ f $(T,a,p)B(a,p)u(a)da, (4) 
Jo 

where $(T, 0, p) is the transition matrix for the system 
X[t,P) = A(t,P)X(t,P). The goal is for the terminal 
and target states to coincide in the function space Hk, so 
setting X(T, j3) = X F ((3), pre-multiplying by $(0, T, p) 
and rearranging yields the integral operator equation 

(Lu)(p)= [ $(0,*,p)B(a,p)u(a)da 
Jo 

= *{O,T,0)X F (fi)-X a (p), (5) 

which implicitly defines the solution it. It is possible to 
use a spectral decomposition, called a singular system, 
of the operator L to find an infinite eigenfunction series 
expansion for the function u G T-Lt of minimum norm 
that satisfies (5). In practice, one can truncate the series 
and approximate the eigenfunctions numerically to esti- 
mate a solution that results in \ \Xt — Xf\\k < e. The 
necessary and sufficient conditions for ensemble control- 
lability are summarized and the method used to numer- 
ically approximate solutions to Fredholm operators is 
presented in detail in Section 4. 

The controllability conditions for general, possibly 
nonlinear, ensemble control problems are presently un- 
known, and analytical control design methods remain a 
challenging problem, although analytical solutions exist 
for certain specific systems. An alternative direction is 
to use a powerful method for numerically solving opti- 
mal control problems called the pseudospectral method. 
This approach has been widely utilized to solve standard 
optimal control problems throughou t the past decade 
for u se in numerous applications (jRoss and FahrooL 
2003), and has recently been exte nded to solve prob- 
lems in optim al ensemble control ()Ruths and Lit 1201 It 
iLi et al.l . [20Tl . The work presented here is focused on 
the derivation of analytical solutions, and numerical ap- 
proximations, to optimal controls for linear stochastic 
ensemble systems. 

To investigate the effect of stochastic factors on the con- 
trollability of ensembles, we extend the notion of en- 
semble controllability to incorporate systems where the 
state is a random variable on a given probability space. 
Consider a modification to the collection (1) by defining 
a family of stochastic control systems 

dX(t, p, u) = A(t, p)X(t, p, w)dt 

+ B(t,f3)u(t)dt + G(t,(3)dS(t,Lj), (6) 

for X G M C E", P G K C K d , u G U C W n , and 
S G V C M fc , where A{t,p) G R" xn , B(t,p) G M" xm , 



and G(t,P) G W ixk have elements that are real L M , 
L 2 , and L 2 functions, respectively, defined on a com- 
pact set D = [0,T] x K C M 2 , and are denoted A G 
L n ^ n {D), B G L£ xm (D), and G G L 2 lxk (D). The equa- 
tion (6) is stated in the Ito sense, where the term dS 
is the differential of a continuous-time stochastic pro- 
cess S : [0, T] x fl — > V C K fe on a probability space 
(fl, T 1 P) . The ensemble controllability conditions for the 
systems (6) depend on the existence of an open-loop con- 
trol u G T-Lt that can steer the instantaneous expected 
state of the stochastic ensemble £ X : D — >• M between 
any points of interest in the function space Hk, where 
£ denotes expectation over the stochastic space f2. Our 
analysis of stochastic ensemble controllability of the col- 
lection (6) is based on the approach taken in the de- 
terministic case, hence a necessary condition is the en- 
semble controllability of the corresponding deterministic 
system, i.e. when G is identically zero. This requirement 
gives rise to the following definition. 

Definition 2. We say that the family of systems (6) is 
ensemble controllable on the function space %k if for 
all e > 0, and all Xo,Xp G Hk, there exist T > 
and an open-loop piecewise-continuous control u G T-Lt 
such that starting from X(0, p) — Xq(P), the final state 
X(T,P) = X T (P) satisfies \\£X T - X F \\ K < e. 

If the system (6) is ensemble controllable, then a natural 
objective for the problem of steering the ensemble from 
Xq to Xp in time T is to minimize the error in the mean 
of the terminal state with respect to the target state, 

Ji = \\SXt — Xf\\k 

= ^JjX T (P)-X F (P)\\ld^p^ ' , (7) 

where || • H2 denotes the Euclidean norm. In addition to 
error in the mean, another important factor in stochastic 
control is MSE, which is minimized using the objective 

J 2 =£\\X T -X F \\\ (8) 
= f £\\X(T,p)-X F \\ 2 2 dv(P), 

J K 

where we invoke Fubini's theorem to evaluate the expec- 
tation before integrating with respect to p. 

Remark 3. If the A, B, G, Xq, and X F all depend con- 
tinuously on the parameter P, then the norm can be 
used in Definition 2 and the objectives J\ and J2 ■ In that 
case, the objective J\ is endowed the physical meaning 
of guaranteeing that the error in the mean is bounded 
irrespective of P G K, and the objective J2 guarantees 
that the MSE is uniformly bounded with respect to p. 
In order to broaden the applicability of our approach, 
we have chosen to use the L2 norm, at the cost of uni- 
formity in this respect. 
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In the following section, we will prove that the optimal 
controls that minimize objectives (7) and (8) coincide, 
and also solve the deterministic control problem when S 
is a Brownian diffusion process. 

3 Stochastic Ensemble Control 

In this section, we derive optimal ensemble controls for 
systems with additive diffusion and jump processes, 
which we model using Brownian motion and Poisson 
counters, respectively. The solution in each case is ob- 
tained implicitly in the form of a Fredholm integral 
equation of the first kind, so it is possible for specific 
cases to obtain an analytical expression of the control 
function as an infinite series of weighted eigenfunctions 
of t he Fre dholm operator as in the deterministic case 
(EL 1201 lh . We present the derivation for the case of 
Brownian diffusion first, and discuss the modifications 
necessary for considering a Poisson type jump process. 

3.1 Brownian Diffusion 

Consider a finite-dimensional time-varying stochastic 
linear system with an additive Brownian noise process 
governed by the Ito differential equation 

dX(t, p, u) = A(t, P)X(t, p, u))dt 

+ B(t,p)u(t)dt + G(t,P)dW(t,u), (9) 



d , u G U C M m , and 
B(t,P) e R" xm , and 



for X G M c M", P G K C 
W G R fe , where A(t,P) G W l> 

G(t,P) G K nxfc have elements that are real Loo, Lo, 
and L2 functions, respectively, defined on a compact set 
D = [0,T] x K C M 2 , and are denoted A G L^ n (D), 
B G L2 Xm (D), and G G L" xfc (L>). The vector valued 
stochastic process W : [0, T] x fl — > R fc consists of in- 
dependent identically distributed (i.i.d.) standard nor- 
mal random variables on the probability space (fi, J 7 , P) , 
and has the natural filtration so that £W(t) = and 
£[W(t)W'(s)] = 7min(t, s). We use £ to denote expec- 
tation over the space SI with respect to measure P, and 
we omit w as an argument of X and W from now on in or- 
der to adher e to a straightfo rward notation for stochas- 
tic calculus ()Brockettl . Il983() . We first derive the control 
that minimizes J\ in (7) as follows. 

Proposition 4. For ensemble controllable system (9), 
the control u G Ht that minimizes J\ , the norm of the 
error in the mean of the terminal state (7), while steering 
the system from the initial state X(0, p) — Xq(P) to the 
target state X(T,P) — Xp{P) must satisfy (5). 

Proof. Using Fubini's theorem and the fact £W(t) = 0, 



hence the expectation of the terminal state is given by 
£X(T,P) = $(T,0,P)X (P) 

+ f <t>(T,a,P)B(<j,P)u(a)da, (11) 



G(t,p)dW(t) 



G(t,P)£dW(t) = 0, (10) 



which coincides with the solution (4) in the deterministic 
case. Setting £X(T,P) = X F (P) results in (5). 



This proposition implies that if (9) is ensemble control- 
lable then one can achieve J\ < e, so that the mean of 
the terminal state approximates the target in the space 
Ti.K- The minimization of the MSE of the terminal state, 
given by objective (8), is addressed in the following the- 
orem, which is the main theoretical result of this paper. 

Theorem 5. For ensemble controllable system (9), 
the control u G %t that steers the system from 
the initial state X(0,/3) — Xq(P) to the termi- 
nal state X(T,P) — Xp(P), while minimizing Ji, 
the mean square error in the terminal state (8), 
must satisfy (5) for all P G K. The minimum 
value of J2 is J K tr C(T, /3)d/i(/3), where C(t,P) — 

f*$(t,o-,P)G(o-,P)G'(o-,p)$'(t,c-,p)do- and $(i,0,/3) 
is the transition matrix for X(t, P) = A(t, P)X(t, P). 

Proof. Let $(t, 0,/3) denote the transition matrix of 
A(t,P). The mean square error in the terminal state 
Xt(P) — X(T, P) as a function of P can be expressed as 

£\ X'y — Xp I2 = £X'rpX r p — 1Xp£X r p + X'pXp , (12) 

where we omit the argument P in order to simplify nota- 
tion, and ' denotes the transpose. Noting that £ X'X = 
tr£XX', where tr denotes the trace of a square matrix, 
we set Z = XX' and apply Ito's rule to obtain 

dZ = (AZ + ZA' + BuX' + Xu'B' + GG')dt 

+ (GX' + XG')dW, (13) 

and taking the expectation results in the deterministic 
linear matrix equation 

E = AT, + HA' + Bu(£X)' + £Xu'B' + GG', (14) 

where £ = £ Z = £XX'. Solving (14) gives rise to 

£(t) = $(t,0)E(0)$'(t,0) 

+ / ^(t,a)B(o-)u(o-)£X'(a)<i>'(t,a)da 
Jo 

+ I ^(t,a)£X(a)u'(a)B'(a)^'(t,a)da 

+ [ $(t,o-)G(o-)G'(o-)$'(t,o-)do- (15) 
Jo 
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as a r esult of elementary linear system theory (Brockctt, 
1970). Observe that the quantities 

C(t)= [ 3>(t,a)G(a)G'(a)$'(t,a)da (16) 
Jo 

and $(t, 0)E(0)$'(t, 0) are independent of the control u. 
Let $t denote $(T, 0), and define a change of variables 

Y(a)= f $(0,r)B(r)M(T)dT, (17) 
Jo 

so that Y(0) = 0. Substituting (11) for £X(cr) in equa- 
tion (15) leads to 

E(T) = $t/ y(a)(x +y(o-))'do-^ 

Jo 

+ $ T / (X + r(cr))y'( C r)d ( T$^ 
Jo 

+ $ T E(0)$y + C(T). (18) 
This, together with F(0) = 0, yields 

£X' T X T =tr [ (X + Y{a)y<P' T <$> T Y{(T)da 
Jo 

+ tr / Y'(<j)& T <S> T (X + Y{a))da 
Jo 

+ tr $ T E(0)$^ +tr C{T) 
= 2X^' T $ T (Y(T) - F(0)) 

+ r'(r)$^$ T y(r) - y'(o)^$ r ^(o) 

+ Xo<f>' T <f> T X + tr C(T) 

= 2x $^$ T r(T) + r'(T)$^$ T r(T) 

+ X $!r$ T Xo + tr C(T). (19) 

Because X is non-random, tr[$TE(0)$T] = H'&T^olll' 
so we can compute the mean square error by substituting 
(19) and (11) into (12) to obtain 

£||-X"t — Xp\\l 

= tr (£(T)) - 2X' F $ T (X + Y{T)) + X' F X F 
= 2X' Q <S>' T <S> T Y{T)+Y'{T)<S> T <S> T Y{T) 

- 2X F $ T (X + Y(T)) + X' F X F 

+ Xfo' T $ T X + tr{C(T)) 
= \\$ T Y(T) + $ T X -X F \\l+trC(T). (20) 

Equation (20) attains its minimum with respect to u 



when Y(T) + X a - ^ X X F = 0, that is, when 

J $(0,T,/3)i?(T,/3)u(T)dT 

= $(0,T,p)X F (p)-X o (p), (21) 

which is the condition (5). When this is satisfied for 
all /3 G K, then the objective J% is minimized at the 
minimum value 

minj 2 = / tr C(T,P)dn(J3). (22) 
Jk 

Remark 6. Theorem 5 together with Proposition 4 
demonstrate that the control minimizing the error in 
the mean of the terminal state also minimizes the mean 
square error in the terminal state for all systems in the 
ensemble (9). The triangle inequality yields — 
Xf\\k < £\\X F ~ £Xt\\k + \\£Xt— Xf\\k, where equal- 
ity is achieved when £X F — X F = in the function space 
U K - In that case £\ \X T ~ X F \ \ K = £\\X T - £X T \\ K . It 
follows that minimizing J7i also minimizes Ji . 



3.2 Poisson Jump Process 

Consider a finite-dimensional time-varying stochastic 
linear system with an additive Poisson jump process 
governed by the Ito differential equation 

dX(t, p, C) = A(t, 0)X{t, P, w)dt + B(t, [3)u(i)dt 

+ G(t,(3)dN(t,0 (23) 

for X G M C R n , j3 G K c R, u G U C K"\ and 
W G M fc , where A(t,j3) G R nXn , B(t,/3) G R nXrn , and 
G(t,/3) G K nxfc have elements that are real ^oo, -^2, 
and L12 functions, respectively, defined on a compact set 
D = [0,T] x K C M 2 , and are denoted A G L 7 ^ n (D), 
B e L% xm (D), and G G L% xk (D). The vector valued 
stochastic process N : [0, T] x EE — > S. k consists of i.i.d. 
Poisson counters on the appropriate probability space 
(S, t?,N), with a vector A £ i' of constant intensities, 
and with the natural filtration so that £N(t) = Xt and 
Cov (N(t),N(s)) = Amin(t,s). Here A denotes the di- 
agonal matrix containing the elements of A in its main 
diagonal, and zero elsewhere. We use £ to denote ex- 
pectation over the space EE with respect to measure N, 
and we omit £ as an argument of X and N from now 
on in order to adhe re to a st r aig htforward notation for 
stochastic calculus (jBrockettl . Il983f) . We first derive the 
control that minimizes J\ in (7) as follows. 



Proposition 7. Given the ensemble system (23), the 
control u G %t which minimizes J\, the error in the 
mean of the terminal state ( 7), while steering the system 
from the initial state X(0, ft) = Xo(/3) to the target state 
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X(T,ft) = X F (fi) must satisfy 



(flu) (/3) = / $(0, a, P)(B(a, P)u(a) + G(a, (3)X)da 
Jo 

= $(0,T,/3)X F (/?)-X (/3). (24) 



Proof. The property £N(t) = Xt and Fubini's theorem 
result in 



G{t,/3)dN(t) 



G(t,f3)Xdt, (25) 



hence the expectation of the terminal state is given by 
£X(T,*3) = $(T,Q,j3)X ((3) 

+ { ^(T,a,(3)B{a,(3)u(a)da 



+ / $(7>,/3)G(cr,/3)AckT. (26) 
Jo 

Setting X(T, (3) = X F (f3), pre-multiplying by $(0, T, 0) 
and rearranging yields integral operator equation (24). 



This proposition implies that if the family (9) is ensemble 
controllable then it is possible to achieve J\ < e, so that 
the mean of the terminal state approximates the target 
state in the space T-Lk ■ The minimization of the mean 
square error of the terminal state, given by objective (8), 
is addressed in the following theorem, which is analogous 
to the main result. 



Theorem 8. Given the ensemble controllable sys- 
tem (23), the control u G T-Lt which steers the 
system from the initial state X(0,/3) = Xq(/3) to 
the terminal state X(T,f3) — X F (f3), while mini- 
mizing Ji, the mean square error in the terminal 
state (8), must satisfy (24) for all f3 G K. The 
minimum value of Ji is J K trC(T, /3)d/i(/3), where 

C(t,(3) = f^(t,a,p)G(a,f3)AG'(a,f3W(t,a,f3)da 
with A = diag(A) being the diagonal matrix of in- 
tensities, and $(£, 0,(3) is the transition matrix for 
X(t,(3)=A(t,p)X(t,f3). 



Proof. We can rewrite (23) in the form dX — AXdt + 
Budt + Yli=i ftdTVi, where gt is the i th column of G, 
and dNi is the i th element of dN. Setting Z = XX' and 



applying Ito's rule for jump processes gives rise to 
dZ = (AZ + ZA + BuX' + Xu'B')dt 

k 

+ J2(9iX' + Xg' i + g i g' i )dN i 

i=l 

= (AZ + ZA + BuX' + Xu'B')dt 
+ GdNX' + XdN'G' + G(diag(dA0)G', (27) 

where diag(diV) G M. kxk is a diagonal matrix with d./V 
on the main diagonal. Taking the expectation results in 
the deterministic linear matrix equation 

S = AY, + TjA + [Bu + GX](£X)' 

+ £X[Bu + GX]' + GAG', (28) 

where S = £XX' and A = diag(A) is the diagonal ma- 
trix of intensities. By applying th e matrix variation of 
constants formula (Brockett. ll970f ). we obtain 

£(£) = $(t,0)S(0)$'(t,0) 

+ f ^(t,a)(B(a)u(a) + G(a)X)£X'{a)^(t,a)da 
Jo 

+ [ <P(t,o-)£X(<T)(B(cT)u(o-)+G(o-)X)'<P'(t,<j)do- 
Jo 

+ [ $(t,o-)G{a)AG'((T)&{t,<T)d<T. (29) 
Jo 

We denote the quantity 



C(t)= $(t,o-)G(a)AG'(o-)$'{t,o-)do- (30) 

and apply a change of variables 

Y(a)= f Q(0,t)[B(t)u(t) + G(r)A]dr, (31) 
Jo 

so that substituting (26) for £X{T,/3) in (29) leads to 
(18), as in the case of Brownian diffusion, and the proof 
continues using the steps in the proof of Theorem 5. It 
can be shown that Ji attains its minimum with respect 
to u when Y(T) + X - <^^}X F = 0, that is, when 



^(O,T,f3)[B(T,0)u(r)+G(r,P)X]dr 



= 9>(0,T,P)X F (fi)-X (p), (32) 

which is (24). When this holds for all (3 G K , the objec- 
tive Ji is minimized at value (22) with G given in (30). 

Remark 9. The proof of Theorem 8 together with 
Proposition 7 demonstrate that the control minimizing 
the error in the mean of the terminal state also mini- 
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mizes the mean squared error in the terminal state for 
all systems in the ensemble (23). 

The above theorems form a theoretical foundation for 
ensemble control of stochastic linear systems. In the fol- 
lowing section, we describe a numerical method for solv- 
ing the equations (5) and (24) that enables the synthesis 
of optimal ensemble controls. 

4 Ensemble Control Synthesis 

We characterize ensemble controllability of a family of 
finite-dimensional time-varying linear systems by cer- 
tain co ndition s on the associated input-to-state opera- 
tor (|L~il . 1201 lT) . The system of interest (1) is considered 
in the Hilbert space setting described in Section 2. The 
Hilbert spaces Ht and Hk with inner products defined 
by (2) and (3), respectively, have well-defined addition 
and scalar multiplication operations, and are thus both 
separable. It follows that the operator defined in equa- 
tion (5) satisfies L G B(Ht,Hk), where B(Ht,Hk) is 
the set of bounded linear operators from Ht to Hk ■ As 
a consequence, L has the adjoint L* satisfying 

(Z,Lu) K = {L*Z,u) T , V£eHif, ueH T - (33) 

We omit the subscript on the inner product when it is 
clear from the context. The conditions characterizing en- 
semble controllability arc related to the singular system 
of the operator L, which requires the following impor- 
tant definition. 

Definition 10. Singular System (|Gohberg et all [20031) : 
Let Y and Z be Hilbert spaces and K : Y — > Z be 
a compact operator. If (ct^,j/„) is an eigensystem of 
KK* and (a^,fi n ) is an eigensystem of K*K, namely, 
KK*v n = a 2 n v n , v n G Z, and K*Kfi n = cj 2 n [i n , fi n € Y, 
where o~ n > 0, and the two systems are related by 

Kfi n = u n v n , K*v n = a n [i n , (34) 

we say that (o~ n , fi n , v n ) is a singular system of K . 

Remark 11. I t can be shown that the operator L in (5) 
is compact (|Li| . l201l[ ). and hence LL* and L*L are both 
compact, self-adjoint, and nonnegative operators. By the 
Spectral theorem, L*L can be represented in terms of its 
positive eigenvalues, namely, L*Ly — J2 n a n(y^ UnjUn 
for all y G Ht- Moreover, because L*L[i n = of^ n , the 
relations L^„ = <J n v n &ndL*v n = a n [i n follow by taking 
v n = (l/a n )Lfi n . This is the infinite dimensional ana- 
logue of the singular value decomposition of a matrix. 

The following theorem provides necessary and suffi- 
cient conditions for ensemble controllability of finite- 
dimensional time- varying linear systems with input-to- 
state operator L, and is valid for stochastic systems 
when the appropriate compact operator is interposed. 



Theorem 12. (fLl . |2011|) The family of systems (1) is 
ensemble controllable on the function space Hk if and 
only if for any given initial and final state Xq Xp G Hk, 
at time t — and t — T < oo respectively, and for 
£ = $(0, T, f3)Xp — Xq, the conditions 

(i) £fe^<oo, (ii) £ett(Z), (35) 

n=l °™ 

hold, where L is the input-to- state operator of the system 
(1) defined in equation (5), with TZ(L) denoting the clo- 
sure of the range space of L, and the collection of triples 
v n ) is a singular system of the linear operator L. 
Moreover, the control law 

oo ^ 

u = V — (£,f„)M„ (36) 

^ On 
n—1 

satisfies (u,u) < (uq,uq) for all uo G IA and uo ^ u, 
where U = {v \ Lv = £ with (i) and (ii)}. 

The above theorem, together with Definition 10, enables 
a method f or approximating solu tions to (5) of mini- 
mum norm (jZlotnik and Lit |2012|) . The integral opera- 
tor equation Lu — £ can be approximated by a linear 
matrix equation, so that the singular system (er„, /i„, u n ) 
of L can be approximated by the singular value de- 
composition of this matrix. Let {/?_/} be a collection of 
points that uniformly distributed throughout the space 
K for j — 0, 1, 2, . . . , P, and let {tk} be a collection of 
points that linearly interpolate the time domain [0, T] for 
k = 0, 1, . . . , N, including endpoints, with tk — tk-i = S. 
Using this grid of nodes, we approximate 

(Lg)((3) = f ${Q,t,f3)B(t,f3)g(t)dt 




$(0,t,f3)B(t,f3)g(t)dt 



f*^8$(0,t h ,p)B(t k ,P)g(t k ) (37) 

k=l 

for each /3 G {fij}- Hence the action of the operator L on 
a function g G Ht can be approximated by the action of 
a block matrix W G M nPxmJV ; with n x m blocks W jk = 
,5$(0,i fc ,ft-)5(^,/3j), on a vector g G M" lN , with N 
blocks g k — g(t k ) of dimension m x 1. If the SVD of this 
matrix is W = UT,V' , and Uj and Vj are columns of U 
and V, respectively, corresponding to the singular value 
Sj, then WW'uj = sjuj and W'Wv k = sjv k . Therefore 
the SVD (sj, Vj, Uj) of the matrix W approximates the 
singular system (o~j, [ij, Vj) of the operator L, where Vj 
and Uj are discretizations of fij and Vj , respectively. Now 

suppose that £ G M. nP is given by i; k = for a 

function £ G Hk- Then the minimum norm solution 
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g* that satisfies Wg = £ is given by g* — W'z where 
WW'z = £ (|Luenbergerl . 1 19691 ). so that applying basic 
properties of the singular value decomposition yields 



mq i,_ 

a =2^— v 3 

3=1 J 



(38) 



The components of the synthesized minimum norm con- 
trol u* = (ul, • - • > UmY are therefore given by 



Ut. 



q Pa 

3=1 



Sfe+m(j-l) 



"Wfc+m(3'-l) ■ 



(39) 



Remark 13. The time and parameter discretizations 
N and P must be chosen such that nP < mN, so that 
the pair (W, £) represents an underdetermined system 
and therefore a minimum norm and not a least squares 
approximation problem. The number q of eigenfunctions 
used in the approximation is limited by q < P. 

5 Examples and Simulations 

We present examples that illustrate the performance of 
our method. We first consider the control of an ensem- 
ble of harmonic oscillators, which are frequently used as 
approximations for periodic phenomena in wide-ranging 
applications from engineering to physics , with parame- 
ter uncertainty as well as additive noise (jBartlett et all 
l2002tlMirrahmil 120041 ). Then, we demonstrate the con- 
trol of an uncertain time- varying stochastic system, as 
well as a three-dim ensional system o f sign ificance to 
quantum transport ()Stefanatos and LH 1201 ll ). 

Example 14. Consider a two-dimensional ensemble 
system with additive Gaussian process dW, given by 

dX(t, u>) = A(t, u)X(t, uo)dt + Bu{t)At + GdW, (40) 



A{t,u) 



with frequency dispersion uj G [—10,10]. We wish to 
steer this harmonic oscillator using a control u(t) £ IR 2 
from initial state Xq = (1,0)' to target state Xp 
:().()'' at terminal time T = 1. The optimal ensemble 
control satisfies (5), which is solved using the method 
in Section 4. In order to avoid numerical conditioning 
errors, we choose q — 5 in (38) such that the largest 
and smallest singular values used satisfy s\/s mq < 10 4 
(Zlot nik and Lil . [2012f ). It is sufficient to sample w at 
21 equidistant values on the interval [—10, 10] to obtain 
this number of eigenfunctions. We use a discretization of 
40001 points over the time interval [0, 1]. It is straight- 
forward to compute the MSE at the terminal state us- 
ing (16) as tr C(T,oj) = TG'G = 0.05, which is invari- 
ant with respect to uj. The optimal control is shown in 



-u 
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' I 


o" 
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Fig. 1. Example 14, where N = 40001 points on [0,T] and 
P = 21 samples of uj 6 [—10, 10] are used, (a) The op- 
timal ensemble control, (b) A stochastic sample trajectory 
for uj — 7 using the Euler-Maruyama method with step 
h = 5 x 10 -4 , and the noise- free path, (c) The expectation 
and (d) MSE of the terminal state for 400 sample paths are 
close to the theoretical values of (0, 0)' and 0.05, respectively. 

Figure 1(a). For the deterministic system, with G = 0, 
it yields a terminal state error of less than 10~ 3 for all 
uj e [—10,10]. The stochastic sys tem is simulated us- 
ing th e Euler-Maruyama method (jKloeden and Plated . 
1999) with sample trajectories shown in Figure 1(b), and 
the statistics at the final state are shown to closely ap- 
proximate the theoretical result in Figure l(c,d). 



Example 15. Consider the system 

dX(t, uj) = A(t, uj)X(t, uj)dt + Bu{t)dt + GdN, (41) 

where A, B, Xq, and Xp are as in Example 14, uj 6 
[-10, 10], G = (0.05,0.05)', T = 1, and diV is a scalar 
Poisson jump process with jump rate A = 20. The opti- 
mal control is shown in Figure 2(a). It is straightforward 
to compute the MSE at the terminal state using (30) 
as tr C(T,uj) = TXG'G = 0.1. The Poisson jump sys- 
tem is simulated using a 4th order Runge-Kutta method 
and pseudo-randomly generated jump times, with sam- 
ple trajectories shown in Figure 2(d), and the statistics 
at the final state are shown to closely approximate the 
theoretical result in Figure 2(b,c). 
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Fig. 2. Example 15, where N = 40001 time nodes on [0, T] 
and P = 21 samples of u> on [—10, 10] are used, (a) The opti- 
mal ensemble control, (b) A stochastic sample trajectory for 
ui = —10. (c) The expectation and (d) the MSE of the ter- 
minal state for 400 sample paths are close to the theoretical 
values of (0,0)' and 0.1, respectively. 



Example 16. Consider a scalar time-varying system 

dx(t, f3) = - sm(0t)x(t, P)dt + u(t)dt + dW, (42) 

with = [—5,5], T = 1, initial state Xq = 1 and 

target state Xp = 0.2. The transition matrix for the sys- 
tem x(t,f3) = — sin(/3£)x(f, ff) is given by $(t, to, ft) — 
exp f to (— sin(/3cr))dcr = exp[i cos((3t) — i cos(/3io)]- The 
optimal control satisfying (5) is shown in Figure 3(a). 
The expected MSE in the terminal state is obtained us- 
ing (16) as C(T, (3) = r Q T exp[| cos(/3t) - § cos(/3i )]dr. 
The expectation of the final state is close to the target as 
shown in Figure 3(b). For terminal time T = 2 the MSE 
for values of /3 £ [—5, 5] closely matches C(2, j3), and for 
(3 = 2 the MSE for T £ [1,3] closely match C(T, 2), as 
shown in Figure 3(c,d). Each statistic is estimated us- 
ing 100 sample paths. A surface plot of C(T,(3) on the 
region (T, (3) £ [1, 3] x [—5, 5] is shown in Figure 3(e). 



Example 17. A quantu m ensemble transport system 
(jStefanatos and Lil . l201ll) with noisy input is given by 



dX = [A(t, lu)X + Bu] dt + GdW, 



(43) 
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Fig. 3. Example 16, where N = 20001 nodes on [0, T] and 
P = 101 samples on j3 £ [—5,5] are used, (a) The optimal 
ensemble control, with q = 9. (b) The mean in the terminal 
state, (c) The MSE of the terminal state for /3 £ [-5,5] 
when T = 2, and for T £ [1,3] when /3 = 2. The statistic 
computed using 400 sample paths with time step h = 10~ 3 
(x) is compared to the theory (line), (d) Surface plot of the 
theoretical MSE C(T,/3) for (T,/3) £ [1,3] x [-5,5]. 
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where dW is a Gaussian process, uj £ [0.8, 1], a = 0.02, 
T = 10, and Xq = (0,0, 1)' and X F = (0,0,0)' are ini- 
tial and target states, respectively. The optimal control 
satisfying (5) is shown in Figure 4(a), and sample tra- 
jectories are given in Figure 4(b). The statistics at the 
final state are computed from 1000 sample paths gen- 
erated u sing a strong 1.5 order sto chastic integration 
method (jKloeden and Platenl . [1999) , and are shown in 
Figure 4(c,d) to closely match the theory. For terminal 
time T = 10, the MSE for all uj £ [0.8, 1] closely matches 
C (10, uj), which is computed numerically by integrating 
(16), and the mean terminal state is also near the target. 

6 Conclusion 

We have derived the optimal open-loop control that 
guides a family of independent, structurally identical, 
finite-dimensional stochastic linear systems with vari- 
ation in system parameters, between initial and target 
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Fig. 4. Example 17, where N = 40001 time nodes on [0, 20] 
and P = 101 samples of u on [0.8, 1] are used, (a) The op- 
timal ensemble control, (b) Sample trajectories for ui = 0.9. 
(c) The norm of the mean in the terminal state, (d) The 
MSE of the terminal state. The statistic computed using 
1000 sample paths with a time step h = 10 -3 . 

states in function space on the parameter domain. En- 
semble control has been extended to stochastic systems 
with additive Brownian noise and Poisson jump pro- 
cesses. It was shown that the same control minimizes 
both the error in the mean and the mean square error in 
the terminal state with respect to the target in function 
space. The optimal control was obtained by solving the 
Fredholm integral equation associated with the system 
dynamics, and was approximated by using the singular 
value decomposition of a matrix that approximates the 
action of the Fredholm operator. We used Monte Carlo 
stochastic integration to evaluate the statistical proper- 
ties of the controlled example systems, and the results 
closely followed our theory with respect to performance 
objectives. Our work has immediate practical applica- 
tions to the control of dynamical systems with additive 
noise and parameter dispersion, which are of interest in 
various areas such as NMR, MRI, and neuroscience, and 
also makes a fundamental contribution to stochastic 
control theory. The novel concepts we explored lead to a 
rich variety of new stochastic control problems, involv- 
ing uncertainty in the system parameters and including 
optimization with respect to the state, control, and time 
horizon, for which standard methods, such as maximum 
principle and stochastic dynamic programming, may 
not be effective to obtain optimal controls. 
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